      module jacobian
!     Jacobian matrix for Newton method 
      use prec
      use bernstein
      implicit none

      contains

      function evaluate_jac(x,j,a,b,n,beta) result(reslt)
!     Element of Jacobian matrix D(i,j) = \partial F_i / \partial beta_j 
!     Usage:
!     D(i,j) = evaluate_jac(x,j,a,b,n,beta)
      implicit none

      real(dp) :: reslt      

      integer, intent(in) :: j,n ! Which basis function [0,n], order of basis functions
      real(dp), intent(in) :: x  ! What collocation point - there are n+1 collocation points
      real(dp), intent(in) :: a,b            ! Interval of calculation
      real(dp), dimension(0:n), intent(in) :: beta ! Beta coefs from previous Newton iteration
      real(dp) :: y

       y = Bernstein_interpolant(beta,n,a,b,x)

      reslt = bernstein_basis_fun_derivative(1,j,n,a,b,x) &
            - 2*bernstein_basis_fun_eval(j,n,a,b,x)*y

      end function evaluate_jac

      function evaluate_fun(x,a,b,n,beta) result(reslt)
!     Evaluates non-linear function written in implicit form F(xi) = 0 at xi for all basis function j=0,n
!     Koristi se:
!     Fun(i,j) = evaluate_fun(x,a,b,n,beta)
      implicit none

      real(dp) :: reslt      

      integer, intent(in) :: n ! order of basis functions
      real(dp), intent(in) :: x  ! What collocation point - there are n+1 collocation points
      real(dp), intent(in) :: a,b            ! Interval of calculation
      real(dp), dimension(0:n), intent(in) :: beta ! Beta coefs from previous Newton iteration
      integer :: k
      real(dp) :: dydx,y


      dydx = 0.0_dp
      do k=0,n,1 ! we use all betas now
      dydx = dydx + beta(k)*bernstein_basis_fun_derivative(1,k,n,a,b,x)
      enddo

!      y = 0.0_dp
!      do k=0,n,1 ! we use all betas now
!      y = y + beta(k)*bernstein_basis_fun_eval(k,n,a,b,x)
!      enddo
       y = Bernstein_interpolant(beta,n,a,b,x)

      reslt = dydx-y*y-x

      end function evaluate_fun

      end module jacobian
